function velocities

%  comparison of numerical phase and group velocity using explicit method on
%     diff(u,x,x)=diff(u,t,t)  for xL < x < xr, 0 < t

% clear all previous variables and plots
clf
clear *

% set parameters
tmax=1.8;
xl=0;
xr=1;

% loop for phase velocity curves
for ic=1:2

	if ic==1
		n=150;
		m=272;
	elseif ic==2
		n=150;
		m=546;
	end;

	% generate the points along the x-axis, x(1)=xL and x(n+2)=xR
	x=linspace(xL,xR,N+2);
	h=x(2)-x(1);

	% generate the points along the t-axis, t(1)=0 and t(m+1)=tmax
	t=linspace(0,tmax,M+1);
	k=t(2)-t(1);
	lamda=k/h;

	% phase velocity
	kb=linspace(0.01,pi/h,100);
	for ik=1:100
		vph(ik)=abs(2*asin(lamda*sin(kb(ik)*h/2))/(k*kb(ik)));
	end;

	subplot(2,1,2)
	if ic==1

		% plot phase velocity
		plot(kb,vph,'-r')
		hold on
	
		% define axes used in plot
		xlabel('wavenumber','fontsize',14,'fontweight','bold')
		ylabel('phase velocity','fontsize',14,'fontweight','bold')
	
		% have matlab use certain plot options (all are optional)
		axis([0 500 0.6 1.02]);
		set(gca,'fontsize',14);
	
	elseif ic==2

		% plot phase velocity
		plot(kb,vph,'-b')
	
	end;
end;

% add in legend
legend([' m = 272 (\lambda = 0.999) '],[' M = 546 (\lambda = 0.498) '],3);
set(findobj(gcf,'tag','legend'),'FontSize',14,'FontWeight','bold');

hold off

% loop for group velocity curves
for ic=1:2

	% set parameters
	if ic==1
		N=150;
		M=272;
	elseif ic==2
		N=150;
		M=546;
	end;

	% generate the points along the x-axis, x(1)=xL and x(N+2)=xR
	x=linspace(xL,xR,N+2);
	h=x(2)-x(1);

	% generate the points along the t-axis, t(1)=0 and t(M+1)=tmax
	t=linspace(0,tmax,M+1);
	k=t(2)-t(1);
	lamda=k/h;

	% group velocity
	kb=linspace(0.01,pi/h,100);
	for ik=1:100
		vg(ik)=abs((explicit(kb(ik)+0.000001,k,h,lamda)-explicit(kb(ik)-0.000001,k,h,lamda))/0.000002);
	end;

	subplot(2,1,1)
	if ic==1

		% plot group velocity
		plot(kb,vg,'-r')
		hold on
	
		% define axes used in plot
		xlabel('Wavenumber','FontSize',14,'FontWeight','bold')
		ylabel('Group Velocity','FontSize',14,'FontWeight','bold')

		% have MATLAB use certain plot options (all are optional)
		axis([0 500 0 1.1]);
		set(gca,'FontSize',14);
	
	elseif ic==2

		% plot group velocity
		plot(kb,vg,'-b')
	
	end;
end;


% explicit
function q=explicit(kb,k,h,lamda)
q=2*asin(lamda*sin(kb*h/2))/k;